Based on the Eshelby's equivalent inclusion method (EIM) and Hill's theorem on discontinuities of elastic fields across the interfaces, a theory for the determination of the stress intensity factors (SIFs) of arbitrarily oriented interacting cracks under non-uniform far-field applied stress (strain) is developed. As shown in this investigation the EIM proposed by Moschovidis and Mura can be extended for treatment of such problems, but their formulations are quite cumbersome and computationally inefficient. An alternative analytical approach is proposed that is computationally more efficient, and unlike the method of Moschovidis and Mura can easily handle complex problems of interacting inhomogeneities and cracks. It is seen that as the interaction between the inhomogeneities becomes stronger, this method yields results that are closer to the solutions reported in the literature than the solutions obtained using the extended EIM of Moschovidis and Mura, which is developed herein. Problems involving combinations of interacting elliptic and penny shape cracks and inhomogeneities are excellent candidates for demonstration of the accuracy and robustness of the present theory, for which the previous EIM produces less accurate results. Due to the limitations imposed on the existing methods, every reported treatment has been tailored for a certain category of problems, and only uniform far-field loadings have been remedied. In contrast, the present theory is more general than the previously reported theories and it encompasses interacting cracks having a variety of geometries subjected to non-uniform far-field applied stress (strain); moreover, it is applicable to modes, I, II, III, and mixed mode fracture.